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Abstract 

Random tilings or packings in the plane are characterized by a size distribution of individual elements 
(domains) and by the statistics of neighbor relations between the domains. Most systems occurring in 
nature or technology have a unimodal distribution of both areas and number of neighbors. Empirically, 
strong correlations between these distributions have been observed and formulated as universal laws. Using 
only the local, correlation-free granocentric model approach with no free parameters, we construct accurate 
analytical descriptions for disc crystallization, size-topology correlations, and Lemaitre's law. 
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Materials with a cellular structure of domains on mesoscopic scales constitute a large variety 
of industrially important materials (e.g. emulsions, foams, polycrystalline metals, ferromagnets) 
as well as the vast majority of living tissues in multicellular organisms [HE]- A description of the 
geometric and topological properties of the domains is crucial for the understanding of the material. 
Most cellular materials exhibit a degree of disorder, so that sizes, shapes, orientations, or neighbor 
relations of domains are determined through statistical distributions functions. 

Of particular interest has been the correlation between the number of neighbors n of a domain, a 
discrete toplogical property, and its area (in 2D) or volume (in 3D), a continuous geometric property. 
The correlation is intuitive - larger domains have more neighbors - but for 2D systems empirical 
laws have been formulated with a claim to some degree of universality. Lewis' law [3, 4J postulates 
a linear dependence of the average area An of ra-sided cells on n for certain systems, while for others 
nonlinear analogs have been observed [3H7] . For many 2D systems with unimodal area distributions, 
a universal correlation is observed between the coefficient of variation of the area ca and that of the 
neighbor distribution c^, cf. [5] and Fig.[l] Fundamental questions have remained unanswered about 
such laws, in particular (i) whether there is a way to understand them analytically, (ii) whether 
they are valid for any tiling of the plane, or dependent on physical characteristics of the system, and 
(iii) whether they can be understood locally by considering a single domain neighborhood (in the 
spirit of a mean-field model), or whether spatial correlations are important. We show here that all 
of these relations follow from a simple, purely local model, that analytical results can be obtained, 
and that empirical data can be classified according to whether these laws are obeyed. 

Very recently, a local model for 3D systems has been developed incorporating a simple connection 
between size and neighbor topology. The granocentric model [SI E] (GM) of Brujic et al. analyzes 
the solid angles taken up by spheres neighboring a central sphere and, given the size distribution of 
spheres, numerically predicts probabilities of numbers of neighbors n and touching neighbors z. To 
self-consistently reproduce the observed n, the model must modify the total available solid angle 
around a sphere (47r) to 47r — 63, where empirically es ~ 0.327r[9]. We extend the GM to the 2D case 
and show that, under certain general approximations, its formalism becomes analytically solvable 
without free parameters. 

We consider 2D polydisperse hard discs with area probability distribution P{A), normalized to 
the mean A = 1. As any packing can be converted to a polygonal tiling by constructing cells around 
the discs (Fig. [2]), this also describes a tiling with a distribution Pp{A). With our main focus on 
such tilings, we will not model touching neighbors, but only the probability P„ of having n general 
neighbors. A central disc of area Ac is then surrounded by n discs of area A (Fig. [2]), where every 
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disc subtends an angle 



= g{A) = 2arcsin + ^/aJA)^ 



(1) 



around the central disc, resulting in the probability distribution of angles fc{4>) = P{A{4)))(dg~^/d(f)). 
Following [9], we can now compute the conditional probability P{n\Ac) and the neighbor probability 
Pn through 



P{n\A,) 



P„ = y" P{n\A,)P{A,)dA,. 



(2) 



Here, Rc,n{(j)) is the probability of the sum of n angles from ([T| reaching (j) and F{(l)max — = 
_^fc{(t^)d4> expresses the probability of the {n + 1)*'' angle exceeding (j)max- The maximum 
angle (pmax available for neighboring discs is again 7^ 27r because of steric effects. But as Euler's 
theorem [TU] demands n = 6 in 2D (for generic three-fold coordinated edges between neighbors), 
'Pmax = 27r + e can be determined analytically: To leading order in ca ^ 1, after fitting 5 (touching) 
neighbors around a disc, there will be a 50% chance of being able to fit a sixth neighbor, i.e., for 
the unmodified (pmax = 27r we obtain n = 5.5. To fit, on average, 0.5 additional neighbors, we need 
to supply the additional average angle subtended by 1/2 disc, i.e., e = n/6 and (pmax = 137r/6. 
Higher-order terms in ca can be obtained, but we will show that this argument is sufficient in the 
entire range of relevant ca- 
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FIG. 1: Dependence of topological disorder c„ on size disorder ca in random cellular networks. The analytical theory (solid line) and 
GM simulations (open diamonds) are in excellent agreement with each other and empirical data (closed symbols: experiments, open 
symbols: simulations). The inset show that RVP and hard-core (he) RVP tilings do not obey this correlation (triangles), but when 
these structures undergo energy minimization in Surface Evolver (arrows), they approach the GM prediction (stars). The dotted line is 
the power-law fit of [5]. 



FIG. 2: Local disc packing in the modified 2D granocentric model: a disc of area is surrounded by discs of area A subtending angles 
4>. The discs serve as templates for polygonal cell construction. 



Eq. ^ can be evaluated numerically for realistic size and angle distributions. In cellular ma- 
terials, these are usually gamma or shifted-gamma distributions [3 IHl [12] ; WeibuU [13] and other 
unimodal distributions have also been used. We employ a gamma distribution fit to /c(0) to eval- 
uate ^ numerically; the results are indicated by the label " GM simulations" in subsequent plots. 
We stress that we have also tried WeibuU and other fits with the same mean and variance, and 
obtained almost indistinguishable results [H]. 

Beyond numerical evaluation, however, we can derive analytical results by instead choosing a 
normal distribution as fit to fc, again preserving the first and second moments. The quality of 
this approximation can be assessed rigorously through an Edgeworth expansion pLS] and confirmed 
numerically. As detailed in the Supplemental Material, the approximation is found to give very 
accurate results independent of the exact distribution function. We stress that our analytics derive 
from the single parameter ca only, with ca ^ 0.5 for the vast majority of systems encountered 
experimentally or theoretically and ca ^ 0.85 for all instances we could find; thus, expansions in 
Ca as a small parameter are possible. Expanding P{A) around A = 1, the coefficient of variation of 
/c(0) is obtained to 0{c\) from ([T]) as 

C4> = ca/ {[g-'rmMi)) ■ (3) 

With Gaussian f^, Rc,n is Gaussian as well, with n-fold mean and variance, while F is an er- 
ror function. Then, P{n\Ac) in (2) is obtained using exp(— (Ax + BY)eii{Cx + D)dx = 
{^/A)eTi{{AD - BC)/VA^TC^) [T6], yielding P(n|/1,) = <l>„+i(c^, - <I>„(ca, A,), where 

$„(c^,AJ = lerf f ^fL^-)-'^— ^ , (4) 
^' '^^ 2 \V2^a^{cA,A,)J ^ ^ 

with 4>{Ac) = 2 arcsin(l/(l + y/Ac)) and, to leading order in ca, the variance = c\s/{{l + s)^(2 + 

s)), with s = y/A^- 

Linearizing the argument of this error function around its root in A^, and again making use 
of a consistent Gaussian approximation with width ca for P{Ac), we obtain a second integral of 
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the same type, yielding, after further expansion for small ca, the explicit prediction for neighbor 
probability P„ = ^n+i(cA) - ^^(ca), with 



where S = sin(137r/12n) = sm{(f)max/'^n) [Ti] . 

We first interpret these results for discs (without constructing polygonal domains). Plotting the 
predicted P„(ca), we observe a pronounced plateau of hexagonal order in the near-monodisperse 



in simulations of hard-disc packings p!7HT9] , where a certain amount of polydispersity [20l [21] or 
bidispersity [22j is necessary to achieve a random structure. The critical ca (terminal polydisper- 
sity), below which crystallization into ordered domains occurs, has been empirically determined to 
be CA,crit ^ 0.1 [T71 [201 In the limit Cyi — )■ 0, ([5| can be further expanded around n = 13/2, 
which yields, to excellent accuracy, a compact formula for the fraction of hexagonally coordinated 
discs. 



where 7^ = 2tt'^/585. Figure [s] shows that the approximations ([5]),([6]) capture both the numerical 
results (GM) and the empirically observed feature of crystallization (Eq. ^ is indistinguishable 
from ([5]) with n = 6). The point of largest curvature of Pe{cA), where the crystallization plateau 
begins, is the critical polydispersity CA,crit = l/V^ ~ 0.106, in very good agreement with empirical 
values of terminal polydispersity. 

The crystallization threshold divides the regime of nearly-monodisperse objects (where additional 
orientational disorder, absent from the model, could maintain randomness, e.g. in monodisperse 
cellular tilings) from that of polydisperse objects (where size disorder overwhelms positional disorder 
[9]). In the latter regime, polygonal cells constructed around the discs should be characterized by 
the same relations as the discs themselves. As an example, the predicted P„ match experimental 
data 1221 [2S] for neighbor probabilities in cellular tissues very well, see inset of Fig. [s] (note that 
topological asymmetry such as P5 7^ P7 is obtained without skewness information from P{A)). In 
the following, we compare the results of the GM and analytical calculations to empirical size-toplogy 
correlation data from cellular systems. In order to translate the disc areas into polygonal areas, we 
simply assume that a disc of area A with n neighbors is inscribed into a regular n-gon (with area 
Ap = {n/Ti) tan(7r/?2)A), and re-normalize the resulting areas. All area-related quantities in Figs. [T] 
and |4] thus refer to polygonal areas. While the assumption of regular polygons misses some of the 





case of small ca (Fig. [3|, so that Pe ~ 1 



for Ca < CA,crit- This effect of crystallization is well known 



Pe = erf (7/0^) , 



(6) 
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FIG. 3: Neighbor probabilities P4 (triangles up), P5 (circles), Pg (diamonds), P7 (squares), Pg (triangles down) of discs in the GM 
simulations. Solid lines are analytical approximations. The plateau of Pe ~ 1 quantifies the range of crystallization near monodispersity 
with CA,crit = 7/\/3 (dashed vertical line). The inset illustrates that the disc P„ arc good predictors for numbers of cellular neighbors 
in random tilings, here for Drosophila wing tissue cells |23l I24| . 



disorder in the system, for cellular systems with domain energy the difference of domain areas from 
those of regular polygons is typically very small (cf. e.g. [26]). 



which converges rapidly (even for c^i — ?■ 1, truncating the series after the third term is very accurate). 
Figure [l] shows the resulting analytical c„(ca), together with numerical results and various empirical 
data points. Also shown is the power-law fit of |5]. We see that, for ca > ca^ctUi the agreement is 
excellent for a large variety of systems. However, simulations of random Voronoi polygon (RVP) 
tilings p7t [25] significantly disagree with the GM theory and appear to lie on a separate line in 
Cn — ca space (inset of Fig. [T]). 

We suggest that this is caused by the RVP domains' lack of a compact shape, which in turn is due 
to the absence of an interfacial energy functional. Almost all physical cellular systems, by contrast, 
exhibit an energy penalty for domain boundaries. We tested this hypothesis by digitizing some of 
the RVP structures from using them as initial conditions for an energy minimization of the 
structure (under constant domain areas) using Surface Evolver [29] (SE) with a uniform interfacial 
energy density. As the inset of Fig. [T] shows, the topological disorder c„ indeed diminishes during 
this process (due to Tl transitions [30]) until the systems approach the observed general c„ — ca 



The Lewis law correlation can be derived within the GM by computing An = J AcP{n\Ac)P{Ac)dAc/ Pr 



The expansion that yields (6) also gives the variance /i2,n of P„, and thus 



Cn = H2nl'^- We find 



oo 




(7) 



correlation. 
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FIG. 4: The dependence of A„ on n is nonlinear in disc and most cellular systems. The analytical prediction for this nonlinear Lewis' 
law agrees well with simulations and experimental data up to n 10. Random Voronoi polyhedra with 28; or without hard cores I27| 
do not follow this prediction (inset). Lewis' data on cucumber tissue 4 also appears to be in this class. 



using the same algorithm as above, including expansions in small ca and around n = 13/2 
Interestingly, the CA-dependences cancel to leading order, giving the c^-independent approximation 



exp 



(8) 



This is in agreement with our numerical results that show An to be insensitive to ca, as well as 
experimental and other empirical data |11|5] that confirm this distinctly nonlinear form of Lewis' law 
(|8|, see Fig. |4| The systems for which the c„ — correlation fails (Fig. [T]) display a linear Lewis' law, 
clearly disagreeing with the GM theory (inset of Fig. |4]) and again suggesting two distinct classes 
of systems distinguishable by compactness of the individual domains. Surface Evolver relaxation 
of the RVP structures cannot restore the nonlinear Lewis' law, because larger areas are not part of 
the area distribution, which is not changed by SE. 
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FIG. 5: Lemaitre's law, a robust correlation between Pg and the toplogical variance /i2,n. The analytical theory from the first three 
terms of (j9jl is given as solid line, confirmed to great accuracy by many empirical data points and GM simulations (symbols). 



A relation between Pq and fi2,n has been reported as Lemaitre's law [TJ EI], with an empirical 
two-part fit suggested by [32], |33]. The present formalism again provides concise formulae: from 
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(|6|, we obtain the leading-order inversion ^/ca'^ a/tt-Pg? which combined with ([t]) yields 



/X2,„(P6) = fy2k - l)erfc {{2k - ^)^P^ 



(9) 



As Pe — ^ 1, this immediately reduces to /i2,n = 1 — -Pe, in agreement with [321 [33]. As Fig. |5] shows, 
(|9| is a uniformly valid approximation that captures both empirical fits in the reported range of 
Pg ^ 0.3 [221 ES]- Even though (joj) relates two topological quantities, it could only be derived 
through the area distribution. Lemaitre's law holds very robustly, being valid for all experimental, 
simulation, and analytical data. In this sense, all of the analytical information about topological 
disorder in these systems is contained in Pg. 

We have shown that a purely local formalism, a generalization of the granocentric model in two 
dimensions, yields accurate and general predictions for size-topology correlations in cellular matter, 
using only the coefficient of variation of the unimodal area distribution. Further generalizations are 
feasible in future work, including neighbor correlations (Aboav-Weaire law) or higher moments of 
the distributions (e.g. to obtain Ck measures [M]). As a leading-order theory, the present formalism 
does surprisingly well and puts a number of empirical relations on a firm footing. 

Acknowledgments: The authors thank J. Brujic and E. Vanden-Eijnden for a frank exchange of 
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